<?php

  # HOGFilter project (c) 2009 PV <1minus1@gmail.com>
  # http://code.google.com/p/hog-filter/
  #
  # A CLI PHP script that filters and pre-processes phylogenetic trees
  # from "HOGENOM" database to... well... make tham better. :)
  #
  # HOGFilter is a free software distributed under the terms of
  # the GNU General Public License http://www.gnu.org/licenses/

  // including required PHP files
  require_once 'config.php';           // configuration parameters
  require_once 'HOGTree.class.php';    // the Tree class
  require_once 'TreeParser.class.php'; // the TreeParser class

  // reading the species names table from file
  say( "Reading the species names table ... " );
  $bacNames = readNamesTable( $config['NamesTableFile'] );
  HOGTree::$bacIDs = array_keys( $bacNames );
  say( "done\n" );

  // reading trees
  say( "Reading and pre-filtering trees.\n" );
  $treeFiles = makeTreeFilesList( $config['TreesDirectory'],
                                  $config['TreesExtension'] );

  // Beginning of Step 1 of Req.Spec.
  $hogTrees = array(); // the main trees container
  $nmin = $config['n']; // minimal acceptable number of bacterial leaves
  $tn = 0; // current tree number (for progress indication only)
  $tc = count( $treeFiles ); // total trees count (for progress indication only)
  foreach ( $treeFiles as $treeName => $fileName )
  {
    ++$tn;
    say( "* Input file '$treeName' ($tn/$tc):\n" );

    // getting a tree from current file
    say( "  Reading... ");
    $newickTree = file_get_contents( $fileName );
      // a text representation of the tree
    if ( $newickTree === false )
      die( "Can't read from file '$fileName'." );
    if ( empty( $newickTree ) )
    {
      say( " The file is empty.\n" );
      continue;
    }
    say( "Parsing... " );
    $tree = parseTree( $newickTree );
    say( "Labels... " );
    $hogTree = new HOGTree( $tree ); // making a tree with parsed labels
    $tree->destroy(); // freeing the memory...
    unset( $tree );
    say( "done\n" );
    // now we've got a fully parsed tree in $hogTree object

    // first we'll check the tree according to Step 1 of Req.Spec.
    $skip = false; // if $skip==true -- the tree will be ignored

    $nL   = $hogTree->numberOfLeaves();
    $nBac = $hogTree->numberOfBacLeaves();
    say( "  Bacterial leaves: $nBac (of $nL)" );
    if ( !$skip )
      if ( !$hogTree->hasNonBacLeaf() ) // no non-bacterial leaves
      {
        say( " = all. " );
        $skip = true;
      }
    if ( !$skip )
      if ( $nBac < $nmin ) // not anough bacterial leaves
      {
        say( " < $nmin. " );
        $skip = true;
      }
    if ( !$skip )
      say( ". " );
    if ( !$skip )
    {
      $hogTrees[$treeName] = $hogTree; // putting the tree in the container
      say( "Added to S1." );
    }
    else
      say( "Skipped." );

    say( "\n" );

  }
  say( "Number of trees in S1 - " . count( $hogTrees ) . ".\n\n" );
  // finished reading trees
  // End of Step 1 of Req.Spec.

  // Beginning of Steps 2 & 3 of Req.Spec.
  say( "Removing trees having p<$config[p]; " .
       "defining root and removing non-bacterial leaves.\n" );
  $tn = 0; // current tree number (for progress indication only)
  $tc = count( $hogTrees ); // total trees count (for progress indication only)
  foreach ( $hogTrees as $hogName => $hogTree )
  {
    ++$tn;
    say( "Tree '$hogName' ($tn/$tc)... " );

    // filtering by 'p' value
    // Step 2 of Req.Spec.
    $hogTree->calcD_rec(); // calculating 'd' values in the whole tree
    $dList = $hogTree->getDList(); // getting the list of 'd' values
    rsort( $dList ); // sorting it in the reverse order
    list( $dm, $dmt ) = $dList[0]; // 'd_max' from Req.Spec
    list( $dpm )      = $dList[1]; // 'd_predmax' from Req.Spec
    $p = sqrt($dm) - sqrt($dpm);   // 'p' from Req.Spec
    say( "p=" . number_format($p,4) . " - " );

    // Step 3 of Req.Spec.
    if ( $p < $config['p'] )
    { // removing the whole tree and skipping to the next one
      $hogTree->destroy(); // freeing the memory...
      unset( $hogTrees[$hogName] );
      say( "removed.\n" );
      continue;
    }
    say( "passed. " );
    // finished filtering by 'p' value

    // choosing which half of the tree to remove
    $lc_l  = $dmt->numberOfLeaves();
      // total leaves count in the 'lower' half of the tree
    $blc_l = $dmt->numberOfBacLeaves();
      // bacterial leaves count in the 'lower' half of the tree
    $lc_u  = $hogTree->numberOfLeaves() - $lc_l;
      // total leaves count in the 'upper' half of the tree
    $blc_u = $hogTree->numberOfBacLeaves() - $blc_l;
      // bacterial leaves count in the 'upper' half of the tree
    $chooseUpper = ( $blc_u/$lc_u > $blc_l/$lc_l );
      // choose the upper (== remove the lower) half of the tree
    $blc = ( $chooseUpper ? $blc_u : $blc_l );
      // bacterial leaves count in the selected half

    // checking minimal count of bacterial leaves
    if ( $blc < $nmin )
    { // removing the whole tree and skipping to the next one
      $hogTree->destroy(); // freeing the memory...
      unset( $hogTrees[$hogName] );
      say( "Bacterial leaves: $blc < $config[n]. Removed.\n" );
      continue;
    }

    // separating two parts from each other
    $parent = $dmt->getParent();
    $parent->unlinkChild( $dmt );
    $bacTree = ( $chooseUpper ? $hogTree : $dmt );
    $nbTree  = ( $chooseUpper ? $dmt : $hogTree );
    $nbTree->destroy(); // destroying the non-bacterial half to free the memory

    // setting the root to the "breakpoint"
    if ( $chooseUpper )
      $bacTree = $parent;
    $bacTree->becomeRoot(); // not a big deal if the root is already there

    // removing non-bacterial branches
    $bacTree->removeNonBacBranches();

    // putting the tree back to the array (for it's root may have changed)
    $hogTrees[$hogName] = $bacTree;
    say( "Added to S2.\n" );
  }
  say( "Number of trees in S2 - " . count( $hogTrees ) . ".\n\n" );
  // End of Steps 2 & 3 of Req.Spec.

  // Beginning of Step 4 of Req.Spec.
  say( "Removing edges with low bootstrap values.\n" );
  say( "(Numbers in brackets are counts of resulting trees for each tree in S2.)\n" );
  $resTrees = array(); // the resulting trees container
  $tn = 0; // current tree number (for progress indication only)
  $tc = count( $hogTrees ); // total trees count (for progress indication only)
  foreach ( $hogTrees as $hogName => $hogTree )
  {
    ++$tn;
    say( "Tree '$hogName' ($tn/$tc)... " );

    $rtc = 0; // resulting trees count (for progress indication only)

    $forrest = array(); // a temporary stack of current tree parts
    array_push( $forrest, $hogTree ); // putting current tree on stack
    while ( count( $forrest ) ) // while stack's not empty
    {
      $curTree = array_pop( $forrest ); // getting a tree from the top of the stack
      if ( $curTree->numberOfLeaves() < $nmin )
        continue; // the tree is too small, skipping it
      $curTree->removeTransitNodes(); // removing the "transit nodes"
      $bsList = $curTree->getBSList(); // getting the list of bootstrap values
      if ( !count( $bsList ) )
        continue; // no more edges have bootstrap values
      sort( $bsList ); // sorting bootstrap values
      list( $mbs, $mbsNode ) = $bsList[0]; // lowest BS value and it's node
      if ( $mbs < $config['b'] ) // if it's below the threshold...
      {
        $tParts = $mbsNode->separate(); // splitting the tre into two parts
        if ( $tParts === false )
          die( "Something bad happend while trying to split the tree. :(" );
        foreach ( $tParts as $tPart )
          array_push( $forrest, $tPart ); // putting both parts on stack
      }
      else // otherwise (all BS values are above the threshold) ...
      {
        $resTrees[] = $curTree; // sending the tree to the results container
        ++ $rtc;
      }
    }

    say( "($rtc) OK\n" );
  }
  // End of Step 4 of Req.Spec.

  // Sending resulting trees to the output file
  $result = "";
  foreach ( $resTrees as $resTree )
    $result .= $resTree->toNewickCompact() . "\n";
  file_put_contents( $config['ResultFileName'], $result );
  // All done.



  // sends the given text to... well stdout, but can send it somewhere else :)
  function say( $text )
  {
    echo $text;
  }

  // reads the species names table from the given file
  function readNamesTable( $filename )
  {
    $bacNamesFile = file( $filename );
    if ( $bacNamesFile === null )
      die( "Can't read names table file '$config[NamesTableFile]'." );
    if ( trim( $bacNamesFile[0] ) !== "CODE\tSPECIES NAME" )
      die( "Unexpected header line in names table file:\n$bacNamesFile[0]." );
    unset( $bacNamesFile[0] );
    $bacNames = array(); // mapping of abbreviated species names to their full names
    foreach ( $bacNamesFile as $row )
      if ( ! empty( $row ) )
      {
        list( $abbr, $name ) = explode( "\t", trim( $row ) );
        $bacNames[$abbr] = $name;
      }
    asort( $bacNames );
    return $bacNames;
  }

  // takes a directory name and file extension and returns
  // the list of matching filenames
  function makeTreeFilesList( $dir, $ext )
  {
    $treeFilesMask = "$dir/*.$ext";
    $treeFilesList = glob( $treeFilesMask );
    if ( empty( $treeFilesList ) )
      die( "No files found matching mask '$treeFilesMask'." );
    $treeFiles = array();
    foreach ( $treeFilesList as $fileName )
    {
      $treeName = preg_replace( "~\.$ext$~", "", basename( $fileName ) );
      $treeFiles[$treeName] = $fileName;
    }
    return $treeFiles;
  }

  // parses the given newick tree and returns a Tree object
  function parseTree( $newickTree )
  {
    $tParser = new TreeParser( $newickTree ); // parsing the tree
    $trees = $tParser->getTrees(); // getting the result
    if ( $tParser->errorOccured() )
      die( "Tree parser error: \"" . $tParser->getErrorMessage() . "\".\n" );
    if ( empty( $trees ) )
      die( "No trees." );

    // destroying trees except the 1st (if there were many)
    $first = true;
    foreach ( $trees as $tree )
      if ( $first )
        $first = false;
      else
        $tree->destroy();

    return $trees[0];
  }

?>
